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Effective  thermal  integration  of  system  components  is  critical  to  the  performance  of  small-scale  (<1 0  kW) 
solid  oxide  fuel  cell  systems.  This  paper  presents  a  steady-state  design  and  simulation  tool  for  a  highly- 
integrated  tubular  SOFC  system.  The  SOFC  is  modeled  using  a  high  fidelity,  one-dimensional  tube  model 
coupled  to  a  three-dimensional  computational  fluid  dynamics  (CFD)  model.  Recuperative  heat  exchange 
between  SOFC  tail-gas  and  inlet  cathode  air  and  reformer  air/fuel  preheat  processes  are  captured  within 
the  CFD  model.  Quasi  one-dimensional  thermal  resistance  models  of  the  tail-gas  combustor  (TGC)  and 
catalytic  partial  oxidation  (CPOx)  complete  the  balance  of  plant  (BoP)  and  SOFC  coupling.  The  simulation 
tool  is  demonstrated  on  a  prototype  66-tube  SOFC  system  with  650  W  of  nominal  gross  power.  Stack 
cooling  predominately  occurs  at  the  external  surface  of  the  tubes  where  radiation  accounts  for  66-92% 
of  heat  transfer.  A  strong  relationship  develops  between  the  power  output  of  a  tube  and  its  view  factor 
to  the  relatively  cold  cylinder  wall  surrounding  the  bundle.  The  bundle  geometry  yields  seven  view 
factor  groupings  which  correspond  to  seven  power  groupings  with  tube  powers  ranging  from  7.6-1 0.8  W. 
Furthermore,  the  low  effectiveness  of  the  co-flow  recuperator  contributes  to  lower  tube  powers  at  the 
bundle  outer  periphery. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Small-scale  SOFCs  are  becoming  increasingly  attractive  over 
conventional  technologies  in  several  applications.  For  portable 
power,  SOFCs  are  predicted  to  have  higher  energy  densities  than 
batteries  for  applications  that  require  long  operational  times  or 
mission  durations.  In  auxiliary  power  applications,  the  high  poten¬ 
tial  efficiency  of  SOFC  technology  outweighs  their  added  cost 
in  comparison  to  conventional  small-scale  diesel  generators.  For 
many  potential  applications,  small-scale  SOFC  systems  become 
most  attractive  when  operating  at  efficiencies  exceeding  40%-LHV. 
High  system  efficiency  is  only  realized  when  the  SOFC  stack  and 
balance  of  plant  (BoP)  are  effectively  integrated.  Strategic  coupling 
of  sinks  and  sources  of  thermal  energy  ensures  thermally  self- 
sustaining  operation  and  reduces  the  need  to  supply  excess  fuel 
to  the  system  for  preheating  flow  streams.  Identification  of  ther¬ 
mal  sources  and  sinks  within  the  system  necessitates  the  need  for 
a  system-level  SOFC  model  that  accounts  for  thermal  interactions 
occurring  between  components. 
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The  cell-stack  in  small-scale  SOFC  systems  is  typically  in  close 
proximity  to  BoP  components  and  depending  on  packaging  and 
cell  geometry,  introduces  a  strong  thermal  interaction  with  adja¬ 
cent  components.  This  is  often  the  situation  for  tubular-based  SOFC 
systems  where  because  of  the  coupling,  a  change  in  SOFC  per¬ 
formance  alters  process  statepoints  throughout  the  system  and 
vice  versa.  Furthermore,  adequately  capturing  the  physics  of  such 
geometries  necessitates  an  understanding  of  the  convective  flow 
field  within  the  bundle.  Reliable  prediction  of  stack  performance 
in  small-scale  tubular  SOFC  systems  therefore  requires  a  model 
that  is  also  inclusive  of  thermofluidic  interactions  between  system 
components. 

System  design  and  simulation  models  found  in  the  extant  liter¬ 
ature  typically  employ  approaches  that  either  (i)  impose  adiabatic 
boundaries  on  all  system  components  [1],  (ii)  employ  thermody¬ 
namic  models  to  predict  the  required  heat  loss  from  components 
based  on  a  given  inlet  and  outlet  state  [2],  or  (iii)  calculate  com¬ 
ponent  heat  losses  without  any  thermal  coupling  [3-5].  These 
approaches  neither  capture  the  thermofluidic  interactions  between 
components  nor  do  they  quantify  the  effect  thermal  coupling  has 
on  system  performance.  When  operating  at  elevated  temperatures 
(500-1000  °C),  some  degree  of  thermal  interaction  amongst  SOFC 
components  occurs  even  for  instances  where  each  component  may 
be  wrapped  in  insulation.  Significant  thermal  interaction  is  typi¬ 
cally  found  in  high  volumetric  power  density  portable  and  mobile 
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SOFC  systems  in  the  1-10  kW  range.  Additionally,  containment  of 
system  thermal  energy  in  such  systems  becomes  more  difficult  as 
the  ratio  of  surface  area-to-volume  increases.  In  such  systems,  com¬ 
pactness  is  a  premium  and  once  robust  operability  is  achieved, 
bulky  insulation  may  be  foregone  in  favor  of  a  more  compact 
system.  Kattke  and  Braun  [6]  have  previously  generated  system 
models  that  incorporate  component  interactions  for  small,  mobile 
planar  SOFC  systems  using  relatively  simple  thermal  resistance 
approaches.  While  such  approaches  can  be  effective  for  systems 
that  contain  primarily  thermal  interactions  (i.e.,  there  is  limited 
dependence  on  the  flow  field  external  to  the  device),  they  are  more 
limited  for  tubular  geometries  that  have  both  high  degrees  of  ther¬ 
mal  and  fluidic  coupling. 

The  objectives  of  the  current  work  are  to  (i)  create  a  steady-state, 
high-fidelity  cell-stack  model  able  to  predict  performance  varia¬ 
tions  amongst  cells,  (ii)  create  a  thermally  coupled  system  model 
around  the  high-fidelity  cell-stack  to  capture  the  thermal  and  flu¬ 
idic  coupling  throughout  the  system,  (iii)  use  simulation  results 
to  suggest  improved  system  designs,  (iv)  and  provide  cell  perfor¬ 
mance  groupings  to  increase  the  accuracy  of  reduced-order  models. 
The  SOFC  system  model  includes  a  tubular  SOFC  stack,  catalytic 
partial  oxidation  (CPOx)  fuel  reformer,  tail-gas  combustor  (TGC), 
recuperator,  and  all  process  flow  conduits.  In  this  paper,  we  first 
present  the  overall  modeling  approach  and  system  geometry  under 
study.  Next  details  on  the  CFD  computational  domain,  user-defined 
function,  individual  component  models,  and  information  exchange 
between  models  is  given.  The  capabilities  of  the  system  simulation 
tool  are  explored  on  a  prototype  66-tube  SOFC  system  design  with 
an  output  power  range  of  500-1000  W  [7]. 

2.  Modeling  approach 

The  level  of  model  fidelity  applied  varies  throughout  the  system 
model.  Electrochemistry,  anode  fluid  dynamics,  and  heat  transfer 
within  the  tube  bundle  is  captured  through  a  previously  developed 
1  -D  electrochemical  tube  model  [8].  Unlike  planar  stacks  where  rel¬ 
atively  uniform  anode  and  cathode  gas  flows  exist  in  small  length 
scale  flow  channels,  tubular  stacks  posses  relatively  large  flow 
areas  external  to  the  tubular  cells  where  spatial  flow  variations 
can  occur.  The  tube  bundle  geometry  under  investigation  is  shown 
in  Fig.  1,  where  the  anode  gas  flow  is  internal  (tube-side)  to  and 
cathode  gas  flow  is  external  (shell-side)  to  the  cells.  A  large  cath¬ 
ode  flow  area  leads  to  spatial  variations  in  fluid  flow,  temperature, 
and  oxidant  concentration.  The  need  to  accurately  capture  spatial 
property  variations  within  the  cathode  due  to  thermofluidic  inter¬ 
actions  is  the  main  motivation  for  incorporating  a  computational 
fluid  dynamics  (CFD)  model  of  the  cathode  gas  volume  as  shown  in 
Fig.  1. 

Each  tubular  cell  in  the  stack  is  modeled  with  a  one-dimensional 
(1-D)  electrochemical  tube  model  that  is  coupled  to  the  three- 
dimensional  (3-D)  CFD  flow  field.  Modeling  of  the  SOFC  stack  alone 
is  not  sufficient  for  the  purposes  of  stack  design  because  of  a 
large  performance  coupling  between  BoP  and  the  stack.  Feedback 
between  BoP  components  and  the  stack  is  required.  Coupling  BoP 
models  to  the  stack  model  provides  the  needed  feedback  for  a  real¬ 
istic  system  modeling  tool.  The  CFD  computational  domain  includes 
the  cathode  (shell-side)  gas  flow,  a  tail-gas  recuperator,  the  CPOx 
fuel/air  preheat  tube  region,  and  stack  endplates.  The  CPOx  and  TGC 
components  and  their  surrounding  geometries  are  modeled  with 
quasi-one  dimensional  thermal  resistance  models  that  are  coupled 
to  both  the  CFD  and  tube  models.  All  models  are  interconnected 
through  passing  of  thermodynamic  states  and  heat  transfer  rates 
at  intersecting  model  boundaries  (see  Fig.  3).  Details  on  system 
geometry  are  given  in  Section  2.1  followed  by  the  details  of  each 
component  model  and  their  coupling  to  create  the  system  modeling 
tool,  Sections  2.2-2.5. 


y 
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Fig.  1.  Top  view  of  66-tube  stack  arrangement.  CFD  grid  surrounds  the  domain  of 
the  tube  model.  Each  cell  is  modeled  independently  with  the  tube  model.  Central 
tube  is  for  fuel/air  preheating. 


2.1.  System  model  geometry 

The  system  design  under  investigation  is  derived  from  a  patent 
application  by  an  SOFC  developer  for  a  portable  power  system  [7]. 
The  patented  stack  is  of  tubular  geometry  with  36  tubes  arranged  in 
a  hexagonal  grid.  The  system,  depicted  in  the  simplified  diagram  of 
Fig.  2,  can  be  generalized  as  a  centrally  located  stack  surrounded  by 
larger  cylindrical  cans  creating  annular  process  flow  channels.  The 
entire  system  is  wrapped  in  insulation  and  operates  as  described 
below. 

First,  oxidant  for  the  cathode  enters  the  system  through  four 
inlet  tubes  at  the  top  of  the  insulation.  Oxidant  is  then  preheated  in 
a  three-fluid  (oxidant,  cathode  gases,  TGC  exhaust)  annular  recu¬ 
perative  heat  exchanger.  Leaving  the  recuperator,  oxidant  turns 
radially  inward  entering  the  cathode.  Oxidant  leaves  the  cathode 
through  concentric  circle  cut-outs  surrounding  each  tube  in  the 
outlet  tube-sheet.  Fuel  entering  the  system  is  first  sent  through  an 
atomizing  spray  nozzle  mixing  with  air  at  the  upper  outside  edge 
of  the  insulation.  The  fuel/air  mixture  is  preheated  as  it  flows  down 
a  centrally  located  tube  before  entering  the  CPOx  reformer.  Leaving 
the  reformer,  reformate  is  distributed  inside  a  fuel  plenum  located 
directly  beneath  the  inlet  tube-sheet.  Reformate  enters  the  anode 
gas  channels  and  flows  up  toward  the  outlet  tube-sheet.  Anode  and 
oxidant  gases  leaving  the  stack  mix  directly  above  the  outlet  tube- 
sheet  and  enter  the  TGC.  TGC  exhaust  gas  flows  radially  outward, 
turns  90°,  and  flows  down  an  annular  channel  in  the  recuperator. 
Exhaust  gas  leaves  the  recuperator  and  is  funneled  radially  inward 
and  is  expelled  from  the  system  through  conduits  at  the  bottom  of 
the  insulation. 

In  order  to  simulate  a  system  of  more  applicable  power  capacity, 
the  36-tube  patent  application  design  was  scaled  to  a  system  with 
a  nominal  power  output  of  650  W.  The  scaled  650  W  system  began 
with  a  66-tube  bundle  geometry  provided  by  the  developer.  The 
cells  are  arranged  in  a  hexagonal  grid  that  is  symmetric  about  the 
x-  andy-axes  as  shown  in  Fig.  1 .  In  the  symmetric  tube  bundle,  there 
are  19  distinct  tubes  as  labeled  in  Fig.  1.  BoP  components  including 
the  recuperator,  CPOx,  and  TGC  are  scaled  from  the  36-tube  sys¬ 
tem  to  fit  around  the  larger  66-tube  stack.  Scaling  the  recuperator 
gas  channels  began  by  attempting  to  maintain  a  constant  Reynolds 
number  between  the  36  and  66-tube  systems.  Equating  Reynolds 
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Cathode  CPOX  Cathode 

Air  Inlet  Air/Fuel  Air  Inlet 


els. 

numbers  simplifies  to  Eq.  (1 )  as  density  and  viscosity  are  constant 
between  the  two  systems.  Examining  the  air  channel  of  the  recuper¬ 
ator,  the  inner  radius  is  set  by  the  outer  radius  of  the  66-tube  bundle 
leading  to  a  larger  hydraulic  diameter,  Dh,  than  in  the  smaller  36- 
tube  system.  The  66-tube  system  also  operates  with  approximately 
nine  times  more  airflow  (i.e.  increased  velocity,  U)  which  requires 
a  very  large  channel  to  maintain  a  constant  Re  #. 


DhU\ 

66  tubes 

DhU\ 

36  tubes 

As  an  alternative  scaling  method,  the  gas  channel  widths  within 
the  recuperator,  CPOx,  and  TGC  of  the  66-tube  system  are  set  equal 
to  the  channel  widths  in  the  36  =  tube  system.  The  result  is  a  nom¬ 
inal  650  W  system  centered  around  a  66-tube  stack  with  a  system 
architecture  equivalent  to  the  36-tube  system,  Fig.  2. 

Details  of  reactor  sizing  in  the  CPOx  and  TGC  along  with 
the  tubular  cell  length,  diameter,  and  thickness  were  extracted 
from  information  provided  by  the  developer.  The  resulting  sys¬ 
tem  dimensions  are  shown  in  Table  1.  Inlet  flow  conditions  to  the 
nominal  650  W  system  were  also  supplied  by  the  developer. 


Table  1 

System  and  tubular  cell  dimensions,  units  in  [cm]. 


Height 

OD 

Ins.  thickness 

CFD  domain 

12.8 

22.6 

3.0 

CPOx 

6.3 

22.6 

2.4 

TGC 

5.8 

21.6 

3.2 

System 

24.9 

22.3 

2.9  (avg) 

SOFC  cells 

12.5  (active) 

1.1 

- 

2.2.  Computational  fluid  dynamics  model 

The  CFD  software  platform  employed  in  the  model  development 
was  ANSYS®  Fluent®.  It  is  also  the  Fluent  software  that  executes  a 
User  Defined  Function  (UDF)  which  integrates  all  component  mod¬ 
els.  In  the  subsequent  sections,  the  Fluent  domain,  settings,  and 
the  UDF  are  discussed. 

2.2.  t.  Domain 

The  domain  of  the  CFD  model  is  illustrated  in  Fig.  2  which 
includes  the  entire  cathode  and  stack  endplates  along  with  the 
majority  of  the  recuperator,  fuel/air  preheat  flow,  and  system  insu¬ 
lation.  The  outer  diameter  (cathode  surface)  of  every  tubular  cell 
represents  the  boundary  between  the  CFD  and  tube  model.  Fluent 
solves  the  mass,  energy,  momentum,  and  species  conservation 
equations  within  the  computational  domain  of  the  CFD  model. 
Owing  to  the  stack  symmetry  (see  Fig.  1 ),  one-quarter  of  the  system 
is  modeled  with  symmetry  boundary  conditions  applied  at  the  x- 
and  y-axes. 

2.2.2.  CFD  model  settings 

Solid  conduction,  convection,  and  radiation  heat  transfer  mech¬ 
anisms  are  modeled  throughout  the  CFD  domain.  Assuming  all 
gas  species  are  transparent  and  non-participating,  a  surface-to- 
surface  radiation  model  is  most  applicable,  but  surface-to-surface  is 
unavailable  in  Fluent  with  symmetry  boundary  conditions.  Alter¬ 
natively,  a  discrete  ordinates  radiation  model  is  utilized.  With 
all  fluids  set  as  transparent,  discrete  ordinates  acts  as  a  surface- 
to-surface  radiation  model.  To  check  the  validity  of  the  discrete 
ordinates  method,  a  complete  66-tube  CFD  model  with  surface-to- 
surface  radiation  was  simulated.  A  comparison  of  numeric  results 
between  surface-to-surface  and  discrete  ordinates  methods  found 
excellent  agreement. 

Diffusion  occurs  within  the  CFD  domain  as  oxygen  is  reduced  at 
the  cathode  electrode  surface  (outer  diameter  of  cells).  The  flux 
of  02  diffusion  to  the  cathode  surface  is  predicted  with  the  1- 
D  tube  model.  Diffusion  from  the  bulk  cathode  gas  is  modeled 
using  Fickian  diffusion  utilizing  temperature  dependent  binary  dif¬ 
fusion  coefficients,  Dy.  Since  the  cathode  is  a  two  species  mixture, 
Do2,m  =  Dq2,n2  and  D02,n2  =  ^n2,o2  where  m  is  the  mixture.  The 
temperature  dependent  polynomial  fit  to  the  binary  diffusion  coef¬ 
ficient,  Do2,n2»  was  calculated  with  Cantera  [9]. 

Do2,n2  =  1-095  X  10“10  •  T2  +  7.993  x  10“8  •  T  -  1.1559  x  10“5 

(2) 

where  T  is  in  Kelvin. 

A  laminar  flow  solver  is  used  throughout  the  CFD  domain  which 
is  appropriate  because  the  largest  Reynolds  number  in  the  flow  field 
is  estimated  at  approximately  1200. 

Piecewise-linear  temperature  dependencies  of  density,  specific 
heat,  thermal  conductivity,  and  viscosity  are  used  for  each  gas 
species.  Thermal  conductivity  and  viscosity  use  a  mass  weighted 
mixing  law  while  the  density  is  calculated  using  the  incompressible 
ideal  gas  law. 

Solid  materials  used  within  the  CFD  domain  include  a  metal, 
insulation,  and  the  SOFC  tube  material.  Metal  components  include 
the  inlet  and  outlet  tube-sheets,  stack  can  wall,  fuel/air  preheat 
tube,  and  the  recuperator  wall  separating  air  and  exhaust  flows. 
All  metals  are  modeled  as  INCONEL®  600  alloy  with  thermal  con¬ 
ductivity  applied  as  a  piece-wise  linear  function  extracted  from 
manufacturer  data  [10].  Thermal  conductivity  values  range  from 
17.3  to  27.5  Wm-1  K-1  in  the  temperature  range  of  interest.  An 
emissivity  of  0.9  is  applied  to  all  metal  surfaces  and  is  representa¬ 
tive  of  an  oxidized  INCONEL®  600  alloy  [10].  The  insulation  around 
the  system  is  modeled  as  fiberboard  with  a  thermal  conductivity 


K.J.  Kattke  et  al.  /  Journal  of  Power  Sources  196  (201 1 )  3790-3802 


3793 


applied  as  a  2nd  order  polynomial  fit  of  manufacturer  data  [11], 
as  shown  below.  An  emissivity  of  0.9  is  assigned  to  all  insulation 
surfaces. 

kfb  =  2.857  x  10“8  T2  -2.743  x  10“6  7  +  2.002  x  10“2  (3) 

A  common  thermal  conductivity,  10.5  Wm_1  K-1,  is  applied  to 
the  tube  solid  in  the  CFD  model  and  1  -D  tube  model.  A  precise  emis¬ 
sivity  of  SOFC  cathodes  is  uncertain;  therefore,  an  emissivity  value 
consistent  with  the  literature  is  used,  £tube  =  0.8  [12,13].  Sensitiv¬ 
ity  analyses  have  revealed  radiation  heat  transfer  to  be  relatively 
insensitive  to  surface  emissivities  because  all  radiation  exchange 
occurs  within  relatively  small  enclosures. 

2.2.3.  User  Defined  Function  (UDF) 

A  UDF  written  in  C  employing  Fluent  built-in  functions  is  used  to 
thermally  integrate  the  CFD  model  to  the  remaining  SOFC  system. 
All  CFD  boundaries  except  the  outer  insulation  surface  represent 
interfaces  with  either  the  1-D  tube  model,  the  CPOx  model,  or  the 
TGC  model.  The  outer  diameter  of  every  tubular  cell  represents  the 
boundary  between  the  CFD  model  and  the  1-D  tube  model.  The  top 
and  bottom  of  the  CFD  domain  represent  the  boundaries  between 
the  TGC  and  CPOX  models,  respectively.  It  is  the  UDF  that  provides 
the  communication  pathways  between  Fluent  and  the  component 
models.  A  schematic  of  the  required  system  model  connections  is 
shown  in  Fig.  3. 

The  UDF  is  called  within  Fluent  at  the  start  of  every  Fluent 
iteration  and  thermally  couples  all  components  at  their  common 
boundaries.  Passing  of  information  from  one  model  to  another  at 
their  intersections  occurs  via  the  writing  and  reading  of  data  files. 
As  an  example,  the  thermodynamic  state  of  fuel/air  leaving  the  CFD 
domain  is  extracted  with  the  UDF  and  written  to  a  data  file  which 


Fig.  3.  Communication  pathways  required  between  all  component  models.  Passing 
of  information  conducted  by  UDF. 


the  CPOx  model  reads  as  the  inlet  to  the  CPOx  reformer.  While  the 
UDF  relays  boundary  conditions  to  the  CFD  model  every  Fluent  iter¬ 
ation,  executing  the  1-D  tube  model,  CPOx  model,  and  TGC  model 
occurs  every  N  Fluent  iterations  to  update  these  boundary  condi¬ 
tions.  The  current  modeling  strategy  is  to  set  N=  5,  but  this  is  an 
input  parameter  that  could  be  optimized  to  reduce  computational 
time. 

2.3.  1-D  electrochemical  tubular  cell  model 

A  previously  developed  1-D  anode-supported  tube  model  is 
employed  to  model  the  electrochemically  active  cell  regions  [8]. 
The  model  incorporates  electrochemistry,  anode  gas  flow,  and  heat 
transfer  within  the  anode  flow  channel  and  tube  solid.  Gas  dif¬ 
fusion  within  the  porous  anode  is  modeled  using  the  Dusty-gas 
model.  Electrochemical  performance  is  based  on  the  Nernst  equa¬ 
tion  with  cathode  and  anode  activation  losses,  concentration  losses, 
and  ohmic  losses.  Axial  conduction  is  assumed  to  occur  within  the 
relatively  thick  anode  structure  only. 

Each  tube  extends  from  the  bottom  of  the  inlet  tube-sheet 
to  the  top  of  the  outlet  tube-sheet.  Electrochemically  active  cell 
area  is  defined  by  the  tube  area  between  the  tube-sheets.  More 
specifically,  the  tubes  are  electrochemically  inactive  within  the 
space  where  the  tubes  protrude  into  the  tube-sheet  because  tubes 
are  being  supported  by  the  tube-sheet.  In  the  case  of  the  inlet 
tube-sheet,  no  oxidant  reaches  the  cathode  surface  and  because  ion 
conduction  is  primarily  normal  to  the  tube  (not  axial)  the  surface  is 
assumed  inactive.  In  contrast,  tubes  are  assumed  electrochemically 
inactive  within  the  outlet  tube-sheet  even  though  oxidant  flow 
reaches  the  cathode  surface.  This  is  a  valid  assumption  because 
the  3.2  mm  thickness  of  the  tube-sheets  is  small  in  comparison  to 
the  125  mm  active  tube  length.  Modeling  of  the  entire  bundle  is 
accomplished  by  simulating  every  tube  in  the  bundle  with  the  1-D 
tube  model. 

2.4.  Coupling  3-D  CFD  domain  to  1-D  tube  model 

The  approach  to  couple  the  3-D  CFD  cathode  to  the  1-D  tube 
model  is  discussed  in  the  following.  The  coupling  process  is  carried 
out  within  the  UDF. 

2.4.1.  Mapping  3-D  mesh  to  1-D  bands 

In  all  cases  considered,  the  CFD  computational  grid  of  the  cath¬ 
ode  at  the  tube  boundaries  is  finer  than  the  1-D  band  discretization 
used  within  the  tube  model  as  illustrated  in  Fig.  4.  The  first  step  in 
the  mapping  process  is  to  number  each  1-D  tube  band  in  order  of 
increasing  axial  location.  A  common  band  grid  is  used  on  every  tube 
within  the  bundle.  Next,  the  axial  locations  of  the  1-D  band  edges 
are  calculated.  Then,  a  loop  begins  over  all  Fluent  control  volumes 
that  border  a  tube  wall.  In  this  loop,  the  centroid  of  the  control  vol¬ 
ume  is  first  calculated.  The  control  volume  is  then  mapped  to  the 


Fig.  4.  Fluent  discretization  at  SOFC  tube  interface  overlaid  on  1-D  tube  model 
bands. 
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tube  band  whose  edges  bound  the  control  volume  centroid  in  the 
axial  direction.  This  procedure  results  in  a  mapping  between  the 
3-D  CFD  grid  and  the  1-D  bands  of  the  tube  model.  This  mapping  is 
used  to  pass  information  between  the  CFD  and  tube  models. 

2.4.2.  Variable  passing 

As  shown  in  Fig.  3,  the  UDF  passes  the  variables  required  for  a 
complete  coupling  of  mass  and  energy  between  the  CFD  and  1-D 
tube  model.  Fluent  provides  the  thermodynamic  state  of  oxidant 
in  the  control  volumes  adjacent  to  the  tube  wall  boundary  as  well 
as  the  temperature  at  the  tube  surface  to  the  tube  model.  The  tube 
model  provides  the  resulting  heat  and  oxygen  flux  occurring  at  the 
tube  surface. 

All  flow  variables  extracted  from  Fluent  are  first  area-averaged 
amongst  all  control  volumes  within  their  assigned  axial  band,  thus 
creating  area-averaged  1-D  profiles  that  the  tube  model  can  inter¬ 
pret.  Every  FLUENT-passed  variable  is  then  sent  to  the  tube  model 
except  for  the  tube  temperature  profile. 

2.4.2. 1.  Tube  temperature  smoothing.  Tube  temperature  profiles 
must  be  smoothed  before  they  are  applied  as  boundary  conditions 
in  the  tube  model.  Due  to  area  averaging,  the  first  derivative  of 
tube  temperature  is  not  smooth,  that  is  step  changes  in  tube  axial 
conduction  result  if  the  area  averaged  tube  temperature  is  applied 
directly  in  the  1-D  tube  model.  Tube  axial  conduction  profiles  are 
smoothed  by  fitting  an  nth  order  polynomial  to  the  area-averaged 
tube  temperature  data.  An  axial  tube  temperature  profile  extracted 
from  the  nth  order  polynomial  fit  is  applied  to  the  tube  model. 

2.4.22.  Interpolating  1-D  tube  in  CFD  model.  All  1-D  tube  model 
variables  passed  to  Fluent  are  applied  using  linear  interpolation. 
As  an  example,  the  value  of  heat  flux  applied  to  a  particular  Fluent 
control  volume  is  determined  with  the  linear  interpolation  scheme 
illustrated  in  Fig.  5  and  given  by  Eqs.  (4)  and  (5), 

*-*1  =  (4) 
Z2-Z1  0-2  -  Q."  1  1 

Qc'  =  '  (0-2  -  Qi )  +  Q!{  (5) 

where  z\ ,  z2  are  the  axial  locations  of  the  tube  model  band  centers. 
zc  is  the  axial  location  of  the  centroid  located  on  the  tube  wall  face 
of  Fluent  control  volume  c.  Q."  and  Q."  are  the  heat  fluxes  calculated 
via  the  tube  model  and  Qf  is  the  heat  flux  applied  to  control  volume 
c  within  the  CFD  model. 

2.4.3.  Boundary  conditions  in  CFD  model 

The  tube  surface  heat  flux  is  applied  directly  as  a  heat  flux  ther¬ 
mal  boundary  condition  at  each  tube  wall  in  the  CFD  model.  Oxygen 


cathode 
gas  flow 


Fig.  5.  Linear  interpolation  of  variable  passing  from  1-D  tube  model  to  3-D  Fluent 
model. 


flux  is  applied  as  a  sink  of  mass  and  energy  to  each  control  volume 
that  borders  the  tube  walls,  where  each  sink  is  calculated  as: 

h*o  face  •  Aface 

Mass  Source  =  — - [=]  kg  m  3  s  1  (6) 

Vcv 

face  '  ^face  '  ho2,face  o 

Energy  Source  = — — - — - [=]Wm  3  (7) 

Vcv 

where  m"2  face  is  the  interpolated  oxygen  mass  flux  rate  applied  at 
a  control  volume  bordering  the  tube  wall  boundary,  Aface  is  the  area 
of  the  tube  wall  face,  VCv  is  the  volume,  and  h02, face  is  the  enthalpy 
of  oxygen  in  said  control  volume. 

2.4.4.  Electro  chemically  inactive  tube  ends 

The  electrochemically  inactive  tube  ends  within  the  tube-sheets 
are  discretized  and  incorporated  within  the  CFD  domain.  The  con¬ 
ductive  heat  flux  at  the  inlet  and  outlet  of  the  electrochemically 
active  tube  length,  calculated  with  the  tube  model,  is  applied  as 
a  heat  flux  boundary  condition  within  Fluent  at  the  interface 
between  the  active  and  inactive  tube  sections.  Anode  gas  flow  in  the 
non-reactive  ends  is  not  discretized  in  Fluent  and  modeled  within 
the  CPOx  and  TGC  models. 

2.5.  Modeling  of  fuel  gas  processing 

Reforming  of  the  liquid  hydrocarbon  fuel  feed  to  the  system 
occurs  within  the  small  CPOx  unit  that  is  integrated  within  the 
stack.  Depleted  anode  fuel  gases  are  completely  oxidized  in  the 
tail-gas  combustor  which  is  located  at  the  end  of  the  tube  bundle 
opposite  of  the  CPOx  unit.  The  TGC  and  CPOx  devices  (see  Fig.  2)  are 
modeled  individually.  Thermal  effects  within  each  model  are  cap¬ 
tured  with  quasi  1-D  thermal  resistance  models.  Lumped  surface 
temperatures  are  assumed  within  the  thermal  resistance  models 
where  lumped  surface  definitions  are  shown  in  Fig.  6.  The  CPOx 
and  TGC  models  share  common  boundaries  with  the  CFD  and  tube 
model.  At  these  model  interfaces,  the  UDF  extracts  and  exchanges 
the  thermodynamic  state,  thus  coupling  the  CPOx  and  TGC  models 
to  the  SOFC  system  model  (see  Fig.  3).  For  example,  the  state  of  pre¬ 
heated  fuel/air  entering  the  CPOx  reformer  is  first  extracted  from 
the  CFD  model  and  used  as  an  input  to  the  CPOx  model.  The  state 
of  reformate  within  the  fuel  plenum  is  extracted  and  used  as  the 
anode  flow  inlet  condition  for  the  1-D  tube  model.  The  CPOx  and 
the  TGC  models  provide  a  critical  coupling  of  the  fuel  reforming 
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Fig.  6.  Lumped  surface  definitions  in  both  the  CPOx  (a)  and  TGC  (b)  thermal  resis¬ 
tance  models.  Conduction  heat  transfer  is  modeled  through  solids  (expect  CPOx  and 
TGC  components)  yielding  an  inner  and  outer  lumped  surface  temperature  on  solids. 
Only  one  surface  node  is  shown  on  solid  surfaces  for  clarity. 
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Fig.  7.  CPOx  model  thermal  resistance  network. 

and  the  combustion  of  unspent  fuel  to  the  CFD  and  tube  models, 
thereby  capturing  the  full  performance  of  the  overall  SOFC  system. 


surface  i.  A  convective  heat  transfer  coefficient,  h,  common  to  all 
surfaces  i,  within  the  CPOx  control  volume  is  used.  In  order  to 
capture  the  presence  of  radiation  heat  transfer  without  the  added 
complexity  of  solving  non-linear  equations,  an  effective  heat  trans¬ 
fer  coefficient  is  used  (Eq.  ( 1 0)).  A  value  of  1 00  W  m-2 1<-1  is  chosen 
to  represent  heat  transfer  throughout  the  CPOx  domain.  While 
there  is  uncertainty  surrounding  this  value,  the  heat  transfer  coef¬ 
ficient  is  seen  as  a  tuning  parameter  which  can  be  varied  in  order 
to  match  model  predictions  to  experimental  data.  Future  publica¬ 
tions  will  also  investigate  the  sensitivity  of  system  parameters  to 
the  chosen  CPOx  heat  transfer  coefficient. 


fa  =  hconv  +  farad  =  100  WlTT2  K"1 


(10) 


Conduction,  Qc0nd,i»  through  solid  regions  i  is  calculated  as: 


Qcond,i  —  C^i, inner  ^i, outer)  " 


M, 


lavg,i 


-1 


(11) 


where  T2iinner  and  Tiou ter  are  the  surface  temperatures  of  solid  region 
i,  Aavg j  is  the  average  surface  area  between  the  inner  and  outer 
surfaces  of  solid  region  i  and  L2  is  the  thickness  of  the  solid  region. 

Reformate  gas  flow  within  the  fuel  plenum  cavity  and  exhaust 
channel  gas  flow  is  assumed  to  be  perfectly  mixed.  Because  gas 
flows  entering  the  CPOx  domain  are  not  at  the  perfectly  mixed 
gas  temperature,  thermal  energy  exchange  is  captured  by  Qg as?2  as 
follows, 


Qgas,i  —  hl/Cpj  Igas ,i) 


(12) 


2.5.1.  CPOx  model 

The  CPOx  reformer  is  a  porous  disc  with  a  catalyst  coating  that 
reforms  the  preheated  fuel/air  mixture  leaving  the  CFD  domain. 
A  gaseous  mixture  of  fuel,  n-hexadecane,  and  air  is  converted  to 
syngas  via  the  overall  reaction  below. 

Ci6H34  +  8  (02  +  3.76N2)  4^  16CO  +  17H2  +  30.0N2  (8) 

It  is  essential  to  couple  a  CPOx  model  to  accurately  predict  sys¬ 
tem  performance.  As  an  example,  conditions  within  the  CFD  model, 
such  as  extent  of  fuel/air  preheat,  directly  effect  the  temperature  of 
reformate  leaving  the  CPOx  which  subsequently  enters  the  anode 
of  the  tube  model.  The  CPOx  model  domain  extends  from  the  bot¬ 
tom  of  the  inlet  tube-sheet  to  the  bottom  of  the  system  insulation 
(see  Fig.  2)  and  also  includes  anode  gas  flow  through  the  non-active 
tube  lengths  within  the  inlet  tube-sheet,  as  discussed  in  Section 
2.4.4.  A  quasi  1-D  thermal  resistance  model  created  in  Python  [14] 
is  applied  to  the  CPOx  domain,  as  shown  in  Fig.  7,  to  capture  ther¬ 
mal  and  fluid  interactions  occurring  in/around  the  CPOx  reformer. 
Reformate  leaving  the  CPOx  reformer  is  assumed  to  be  in  chemical 
equilibrium  [15,16].  Within  the  thermal  resistance  model  is  a  CPOx 
reformer,  fuel  plenum  wall,  fuel  plenum  gases,  an  inner  can  that 
separates  the  air  channel  from  the  exhaust  channel,  air  and  exhaust 
gases  located  within  the  channels,  and  insulation  as  shown  in  Fig.  6. 

2.5. 1. 1.  Thermal  resistance  model.  Convection  and  conduction  heat 
transfer  are  modeled  within  the  CPOx  region  assuming  lumped 
surface  temperatures  on  all  components  (see  Fig.  6).  The  surface 
temperature  of  the  CPOx  reformer  is  taken  as  the  average  of  the 
inlet  and  outlet  gas  temperatures.  Thermal  coupling  within  the 
CPOx  domain  is  accomplished  through  a  series  of  energy  balances 
on  all  surfaces  and  the  gases  within  the  fuel  plenum  and  exhaust 
channel.  The  convective  heat  transfer,  Qc0nv,;>  from  any  surface  i  is 
calculated  as: 

Qconv,i  =  (j'i  ~  7gas,j)  (9) 

where  T2  is  the  lumped  temperature  and  A2  is  the  area  of  surface 
i  and  Tgasj  is  the  free  stream  temperature  of  gas  in  contact  with 


where  Qgas,2  is  the  amount  of  thermal  energy  added  at  each  gas 
volume  node  i,  (Tl  in  ^gas.i)  is  the  temperature  difference  between 
the  gas  entering  the  volume  and  the  perfectly  mixed  temperature, 
m2  is  the  flow  rate,  and  Cp2  is  the  specific  heat  calculated  at  the 
average  of  the  inlet  and  perfectly  mixed  temperatures.  Anode  gas 
flow  inside  the  inactive  tube  ends  is  added  to  the  fuel  plenum  cavity 
volume.  During  the  physical  steady-state  operation  of  this  system, 
all  air  flowing  down  the  recuperator  enters  the  cathode  and  no  air 
flow  enters  the  air  channel  surrounding  the  fuel  plenum;  therefore, 
there  is  no  thermal  energy  source,  Qgas?2,  in  the  air  channel,  Fig.  7. 

2.5.12.  Model  integration  to  system  model.  The  thermodynamic 
state  of  all  flow  inlets  to  the  CPOx  model  which  include  the  fuel/air 
preheat  and  exhaust  leaving  the  recuperator  are  extracted  from 
Fluent  via  the  UDF.  The  equilibrium  CPOx  reformate  is  sent  as  the 
anode  inlet  condition  for  the  stack.  A  perfectly  mixed  condition  is 
applied  within  the  fuel  plenum;  therefore,  reformate  is  uniform  in 
temperature,  pressure,  and  composition  entering  the  anode  of  all 
tubes  within  the  bundle.  Besides  flow  interfaces,  there  are  also  solid 
interfaces  at  the  CPOx  boundary  between  the  CFD  and  tube  mod¬ 
els.  Adiabatic  boundary  conditions  are  applied  at  the  interfaces  of 
the  system  insulation  and  the  inner  can  wall.  The  top  of  the  fuel 
plenum  cavity  as  well  as  anode  gas  channels  in  inactive  tube  sec¬ 
tions  are  bound  by  surfaces  within  the  CFD  model.  At  these  CFD 
surfaces  a  convective  thermal  boundary  condition  is  applied  with 
h  =  100  Wm-2  K-1  and  a  free  stream  temperature  equal  to  the  fuel 
plenum  gas  temperature  as  calculated  in  the  CPOx  thermal  resis¬ 
tance  model.  The  total  heat  transfer  at  the  interface  calculated  by 
the  CFD  model,  Qc fd,cpox»  is  added  to  the  energy  balance  at  the  fuel 
plenum  gas  node  within  the  CPOx  thermal  model. 

2.5.2.  TGC  model 

Thermal  interactions  in/around  the  TGC  are  captured  with  a 
thermal  resistance  network  following  the  same  principles  used  in 
the  CPOx  model.  Unlike  the  CPOx  model,  the  TGC  model  assumes 
complete  combustion  [  1 ,3  ]  of  any  unspent  fuel  leaving  the  stack.  In 
addition,  the  TGC  must  capture  the  preheating  of  two  process  gas 
streams,  (i)  air  bound  for  the  cathode  and  (ii)  a  mixture  of  fuel  and 
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Table  2 

Simulation  parameters. 


Fuel/air  inlet 

Fuel  type 

C16H34 

r/°c 

40 

P/kPa 

92.65 

Air  inlet 

r/°c 

20 

P/kPa 

84.37 

Stack 

javg/Acm-2 

0.349 

System 

CPOx:  O/C 

1.1 

A.  air 

2.55 

air  bound  for  the  CPOx.  Details  of  the  TGC  model  development  can 
be  found  in  Appendix  A. 

3.  Simulation 

The  objective  of  the  simulation  is  to  demonstrate  the  capabili¬ 
ties  and  utility  of  the  developed,  high-fidelity  SOFC  modeling  tool. 
To  begin,  an  energy  and  mass  balance  are  calculated  to  check  the 
integrity  of  the  system  model.  System  level  results  are  presented 
illustrating  the  capability  of  the  model  to  highlight  ineffective  sys¬ 
tem  architectures.  Next,  a  detailed  stack  analysis  illustrates  the 
non-uniform  performance  within  the  stack  and  the  usefulness  of 
the  model  to  stack  developers. 

The  nominal  gross  electric  power  output  of  the  SOFC  system 
is  650  W.  Thus,  the  electrical  power  target  for  the  simulation  was 
±10%  of  the  nominal  value.  To  vary  bundle  power  without  chang¬ 
ing  inlet  conditions,  the  insulation  thickness  around  each  modeling 
domain  (CFD,  CPOx,  and  TGC)  was  varied  by  a  common  factor.  Uti¬ 
lizing  a  2X  insulation  factor,  the  system  model  resulted  in  637  W 
gross  electrical  power  with  insulation  30,  24,  and  31.5  mm  thick 
surrounding  the  CFD,  CPOx,  and  TGC  domains,  respectively. 

3. t.  Simulation  input  parameters 


Fig.  8.  System  process  statepoint  locations. 


3.2.  Verification  of  system  model  integrity  -  energy  and  mass 
residuals 


The  symmetric  tube  bundle,  shown  in  Fig.  1,  was  simulated 
yielding  1 9  independent  tube  simulations.  Inlet  conditions  and  sys¬ 
tem  parameters  applicable  to  a  650  W  stack  were  supplied  by  the 
SOFC  developer  and  are  summarized  in  Table  2.  The  amount  of  sto¬ 
ichiometric  air,  Aair,  for  this  hexadecane  fueled  system  is  calculated 
using  the  following, 


1  .  _  (  h°2 
alr  \  24.5  ■  nC|6H34 


system  inlet 


(13) 


The  integrity  of  the  system  model  was  verified  by  performing 
mass  and  energy  balances  around  the  system  and  component  mod¬ 
els.  The  resulting  residuals  are  shown  in  Table  3.  The  total  energy 
residual  on  the  system  is  0.50%  of  total  energy  input  into  the  sys¬ 
tem  and  the  combined  CFD  plus  tube  bundle  energy  residual  is  less 
than  2%  of  SOFC  stack  power.  Mass  balances  across  the  TGC  and 
CPOx  models  with  a  slight  imbalance  seen  in  the  CFD  plus  tube 
bundle  model.  This  mass  imbalance  is  very  low  being  less  than  5% 
of  the  relatively  small  02  mass  consumed  in  the  stack. 


Ambient  conditions  of  Tamb  =  20  °C  and  Pamb  =  83  kPa  surround  the 
system  insulation.  The  thermal  boundary  imposed  to  the  outer 
insulation  surface  surrounding  the  system  is  convective  with 
h  =  10W  m-2  K-1 .  In  addition  to  convection,  a  radiation  pathway  is 
imposed  at  the  insulation  periphery  in  the  CFD  model,  but  the  addi¬ 
tion  of  radiation  has  little  to  no  effect  because  of  the  low  insulation 
skin  temperature,  50  °C. 


3.3.  Simulation  of  the  SOFC  system 

The  SOFC  power  generator  is  simulated  operating  on  a  liq¬ 
uid  hexadecane  fuel.  Simulation  results  are  presented  for  both 
the  overall  system  and  for  the  tube  bundle.  Table  4  presents  the 
thermodynamic  states  throughout  the  system  with  statepoints  cor¬ 
responding  to  Fig.  8.  The  energy  content  at  each  statepoint  is  the 


Table  3 

Energy  and  mass  residuals. 


Model 

Energy  residual 

Mass  residual 

Value/W 

%  Total  inlet  energy 

%  SOFC  power 

Value/gs_1 

%  Total  inlet  flow 

%  O2  consumed  in  stack 

System 

14.9 

0.5% 

2.3% 

-3.4E-4 

0.01% 

4.2% 

CFD  +  tube  bundle 

12.4 

0.4% 

1.9% 

-3.4E-4 

0.01% 

4.2% 

TGC  region 

1.9 

0.1% 

0.3% 

0.0 

0.00% 

0.0% 

CPOx  region 

0.6 

0.0% 

0.1% 

0.0 

0.00% 

0.0% 
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Table  4 

System  statepoints. 


Statepoint  T/°C  P/kPa  Flowrate/gs_1  Molar  composition  E/W 


h2 

CO 

ch4 

C16H34 

co2 

h2o 

o2 

n2 

1 

20 

84.37 

2.410 

0.0 

0.0 

0.0 

0.00 

0.0 

0.0 

0.21 

0.79 

0 

2 

111 

84.37 

2.410 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.21 

0.79 

223 

3 

668 

84.25 

2.410 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.21 

0.79 

1664 

4 

699 

84.22 

2.330 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.19 

0.81 

1695 

5 

40 

92.65 

0.402 

0.0 

0.0 

0.0 

0.02 

0.0 

0.0 

0.21 

0.77 

3010 

6 

294 

92.65 

0.402 

0.0 

0.0 

0.0 

0.02 

0.0 

0.0 

0.21 

0.77 

3157 

7 

422 

92.65 

0.402 

0.0 

0.0 

0.0 

0.02 

0.0 

0.0 

0.21 

0.77 

3227 

8 

1218 

92.65 

0.402 

0.24 

0.24 

0.0 

0.0 

0.01 

0.02 

0.0 

0.50 

3226 

9 

816 

84.22 

0.402 

0.24 

0.24 

0.0 

0.0 

0.01 

0.02 

0.0 

0.50 

2974 

10 

718 

84.22 

0.483 

0.10 

0.10 

0.0 

0.00 

0.14 

0.16 

0.0 

0.50 

1633 

11 

716 

84.22 

2.813 

0.02 

0.02 

0.0 

0.0 

0.03 

0.03 

0.15 

0.76 

3371 

12 

999 

84.22 

2.813 

0.0 

0.0 

0.0 

0.0 

0.05 

0.05 

0.13 

0.77 

3345 

13 

868 

83.27 

2.813 

0.0 

0.0 

0.0 

0.0 

0.05 

0.05 

0.13 

0.77 

2889 

14 

656 

83.22 

2.813 

0.0 

0.0 

0.0 

0.0 

0.05 

0.05 

0.13 

0.77 

2177 

15 

682 

83.00 

2.813 

0.0 

0.0 

0.0 

0.0 

0.05 

0.05 

0.13 

0.77 

2264 

Power 

637 

Heat  loss 

94 

summation  of  (i)  thermal  energy  released  when  cooled  to  the  ambi¬ 
ent  temperature,  (ii)  chemical  energy  released  with  the  oxidation  of 
any  fuels  present,  and  (iii)  latent  energy  associated  with  condensing 
water  if  present. 

System  air  is  preheated  from  20  °C  to  668  °C  before  entering  the 
cathode  at  statepoint  3.  The  fuel/air  gas  mixture  is  preheated  from 
40  °C  to  422  °C  (statepoint  7)  prior  to  entering  the  CPOx  reformer 
and  subsequently  enters  the  anode  at  816  °C  (statepoint  9).  The 
bundle  at  an  average  temperature  of  745  °C  produces  637  W  of 
gross  power  at  43.3  V  and  14.72  A.  TGC  exhaust  gases  enter  the 
recuperator  at  868  °C  (statepoint  13),  leaving  the  recuperator  at 
656  °C  (statepoint  14),  and  finally  are  heated  slightly  while  flowing 
through  the  CPOx  domain  leaving  the  system  at  682  °C  (statepoint 
15).  The  co-flow  recuperator  leads  to  a  high  system  exhaust  tem¬ 
perature.  Of  the  3.01  kW  of  energy  entering  the  system,  75%  is 
convected  away  through  exhaust  gases  with  21%  converted  to  dc 
electrical  power  in  the  SOFC  tube  bundle.  Conductive  heat  loss 
through  system  insulation  only  accounts  for  3%  of  total  system 
energy  input.  Table  5  highlights  the  predicted  operating  conditions 
of  the  system  utilizing  the  following  definitions, 


Up  _  \  H2  )  stack  consumed 

(4ftCH4+ftH2+»C0)an0deinlet 

(14) 

U  _  (  ^2  )  stack  consumed 
(^02  )  cathode  inlet 

(15) 

^DC,  stack 

^yStem“(nfuel'HHVfuel)system  inlet 

(16) 

_  ^DC,  stack 

(^fuel  '  HHVfuei)anocjejnjet 

(17) 

where  UF  is  the  fuel  utilization  and  Uox  is  the  oxygen  utilization 
within  the  stack,  ^system  and  ij Sofc  are  the  system  and  stack  effi¬ 
ciencies,  respectively.  Pdc, stack  is  the  DC  power  output  from  the 


Table  5 

System  operating  conditions. 


Operating  conditions 

^cathode,  avg/kPa 

84.24 

Panode,  avg/kPa 

84.22 

UF 

0.56 

Uox 

0.14 

^SOFC 

25.1% 

system 

21.1% 

stack  and  found  in  both  efficiency  definitions  because  blowers, 
pumps,  and  power  conditioning  are  not  modeled  in  this  study. 
The  relatively  low  system  efficiency  registered  in  this  simulation 
is  consistent  with  SOFC  developer  performance  and  substantially 
affected  by  the  low  fuel  utilization. 

The  predicted  CPOx  reformer  exhaust  temperature,  1218  °C, 
could  lead  to  sintering  within  the  reformer.  Heat  loss  from  the 
reformer  is  directly  coupled  to  the  estimated  reformer  skin  tem¬ 
perature  used  within  the  thermal  resistance  network.  Averaging 
inlet  and  outlet  gas  flow  temperatures  to  estimate  the  CPOx  skin 
temperature  adds  uncertainty  to  the  thermal  model.  Experimental 
temperature  profile  data  of  the  reformer  would  remove  the  uncer¬ 
tainty  surrounding  the  CPOx  skin  temperature  and  predicted  CPOx 
exhaust  temperature. 

3.4.  Tube  bundle  results 

Tube  bundle  performance  is  investigated  with  the  aim  of  under¬ 
standing  tube  temperature  and  oxygen  distributions,  flowfield 
characteristics,  and  identification  of  potential  model  reduction 
methods  by  strategic  groupings  of  tubes. 

3.4 A.  Tube  performance  groupings 

For  every  tube  in  the  bundle,  radiation  is  observed  to  be  the 
dominant  heat  transfer  mechanism  for  stack  cooling  at  the  outer 
tube  surface.  The  annular  stack  can  surrounding  the  bundle  sepa¬ 
rates  the  cathode  gases  from  the  relatively  cold  air  being  preheated 
in  the  recuperator.  The  stack  can  has  a  non-linear  temperature  dis¬ 
tribution  with  a  maximum  of  691  °C,  a  minimum  of  523  °C,  and  a 
621  °C  average  temperature.  With  the  average  bundle  temperature 
at  745  °C,  a  large  temperature  driving  force  for  radiation  heat  trans¬ 
fer  to  the  stack  can  exists.  Tubes  with  larger  radiation  view  factors 
to  the  stack  can  are  observed  to  transfer  a  proportionally  larger 
amount  of  radiation  thermal  energy  causing  their  temperatures  to 
decrease  along  with  tube  power.  Power  disparities  are  observed 
within  the  bundle  as  cell  power  ranges  from  7.6  to  10.8  W  with  an 
average  of  9.7  W. 

Power  disparities  are  largely  due  to  a  strong  functional  relation¬ 
ship  between  the  power  output  of  a  tube  and  the  view  factor  from 
the  tube  to  the  stack  can.  This  relationship  is  shown  in  Fig.  9  where 
the  tubes  are  arranged  in  increasing  power  output  on  the  x-axis. 
Tubes  with  similar  view  factors  to  the  stack  also  have  very  simi¬ 
lar  power  outputs.  For  the  bundle  configuration  under  study,  seven 
view  factor  groupings  exist  which  lead  to  seven  tube  power  group¬ 
ings.  Groupings  6  and  7  could  be  combined  into  a  single  power 
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Fig.  9.  Plot  of  tube  power  and  view  factor  from  tube  to  stack  can  surface  for  every 
tube  in  bundle.  Arranged  in  groupings  of  similar  view  factors. 

group,  but  they  are  left  as  distinct  groups  because  of  differences 
in  axial  temperature  profiles  which  is  discussed  further  in  Section 
3.4.2.  Variations  in  tube  performance  point  to  the  potential  for  sub¬ 
stantial  inaccuracy  in  stack  power  prediction  if  the  performance  of  a 
single  tube  is  extrapolated  to  emulate  the  performance  of  an  entire 
tube  bundle.  Resulting  power  groupings  suggest  at  least  six  if  not 
seven  tube  simulations  are  required  in  order  to  predict  the  perfor¬ 
mance  of  the  stack,  where  each  tube  simulation  requires  a  unique 
set  of  thermal  boundary  conditions. 

The  dominance  of  radiation  in  the  bundle  is  explicitly  illustrated 
in  Fig.  10  where  individual  tube  heat  losses  have  been  averaged 
within  power  groupings.  For  tubes  located  closer  to  the  stack  radial 
center,  outer  radial  periphery  tubes  act  as  radiation  shields  effec¬ 
tively  blocking  the  view  to  the  relatively  cold  stack  can.  This  is 
convincingly  seen  in  Fig.  10  where  the  outer  radial  tubes  (groups 
1  and  2)  have  the  greatest  percentage  of  radiation  loss,  87-92%, 
compared  to  inner  radial  tubes  (groups  6  and  7)  where  66-67% 
of  heat  transfer  is  due  to  radiation.  Because  of  this  shielding,  tube 
power  groupings  are  also  a  function  of  radial  tube  location  as  clearly 
shown  in  Fig.  9  where  tubes  at  the  inner  radii  produce  the  great¬ 
est  power  and  tubes  at  the  outer  bundle  radii  produce  the  lowest 
power. 

Because  the  stack  is  wired  in  electrical  series,  cell  voltage  varies 
throughout  the  bundle  as  shown  in  Fig.  1 1 .  Cell  voltages  vary  from 
0.52  to  0.73  V  at  the  lowest  power  and  highest  power  tubes,  respec¬ 
tively.  Voltages  lower  than  0.6  V  occurring  in  tube  groupings  1  and 


Fig.  10.  Heat  transfer  pathways  in  tube  groupings.  Percentage  of  radiation  written 
on  graph. 


Fig.  11.  Operating  voltage  of  every  tube  within  bundle. 

2  are  of  concern  as  oxygen  ions  may  begin  to  oxidize  Ni  in  the 
anode.  Further  research  needs  to  determine  the  voltage  limit  to 
avoid  Ni  oxidization.  If  tube  groups  are  below  the  voltage  limit, 
each  cell  within  a  grouping  could  be  wired  in  series  and  each 
grouping  could  be  wired  in  parallel.  Thus,  the  current  load  for  each 
grouping  could  be  varied  in  order  to  maintain  acceptably  high  cell 
voltages. 

3.4.2.  Temperature  and  02  axial  profiles 

Area-averaged  axial  tube  temperature  profiles,  T(z),  along  with 
surface  oxygen  molar  concentrations,  X0 2(z),  are  shown  for  every 
tube  in  Fig.  12.  Fig.  12  is  organized  utilizing  the  power  groups 
defined  in  Fig.  9.  As  tube  power  is  a  strong  function  of  tempera¬ 
ture,  tube  temperatures  within  a  given  power  group  remain  within 
1 0  °C  of  each  other.  Conduction  from  the  relatively  hot  CPOx  region 
(81 6  °C  anode  gas  inlet)  combined  with  high  localized  current  den¬ 
sity  contribute  to  the  maximum  cell  temperature  at  the  anode  inlet. 
Because  air  enters  the  cathode  at  a  relatively  cold  temperature 
of  668  °C,  all  tube  temperatures  initially  decrease  with  increasing 
axial  distance.  Due  to  their  close  proximity  to  the  stack  can,  tubes 
in  groups  1  and  2  continue  to  drop  in  temperature  with  a  small 
increase  near  the  anode  outlet  due  to  the  heat  supplied  by  the  TGC. 
All  remaining  tubes  in  the  bundle  produce  enough  thermal  energy 
due  to  irreversibilities  in  the  electrochemical  reactions  to  increase 
their  temperatures  after  the  initial  temperature  decrease.  The  AT 
between  the  bundle  and  stack  can  increases  in  the  direction  of  cath¬ 
ode  and  anode  flow  because  cold  oxidant  enters  the  recuperator 
near  the  anode  and  cathode  outlets.  Radiation  exchange  overcomes 
internal  tube  heat  generation  causing  temperatures  to  decrease  at 
the  local  temperature  maximum  occurring  approximately  halfway 
down  the  length  of  the  bundle. 

Surface  oxygen  concentrations  are  dependent  on  bulk  convec¬ 
tive  cathode  mixing  as  well  as  diffusion.  Oxygen  diffusion  from  the 
bulk  cathode  to  the  tube  surface  is  driven  by  the  current  produced 
by  a  tube.  With  the  total  current,  14.72  A,  of  each  tube  held  con¬ 
stant,  the  sum  of  02  diffusion  to  the  tube  surface  is  constant  for 
each  tube.  Localized  current  densities  are  directly  proportional  to 
tube  temperature  as  seen  in  Fig.  12b.  As  the  temperature  of  group 
4  increases,  02  concentration  decreases  at  an  increased  rate  due 
to  a  higher  localized  current  density.  Interestingly,  surface  02  con¬ 
centrations  do  not  decrease  in  the  flow  direction  along  the  length 
of  the  cell.  The  3-D  cathode  model  captures  increases  in  axial  02 
concentrations  due  to  variations  in  convective  mixing.  For  exam¬ 
ple,  02  concentrations  increase  near  the  cathode  outlet  because 
cathode  gases  exit  through  concentric  circle  cut-outs  surrounding 
each  cell  in  the  outlet  tube-sheet.  Combining  the  quasi-radial  bun- 
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Fig.  12.  Temperature  and  surface  oxygen  molar  concentration  profiles  of  every  tube  in  bundle.  Tubes  grouped  based  on  power  groupings  in  Fig.  9.  (a)  Tube  groupings  1  and 
2.  (b)  Tube  groupings  3,  4,  and  5.  (c)  Tube  groupings  6  and  7. 


die  symmetry  and  the  radial  cathode  inlet  to  the  bundle,  variations 
in  convective  mixing  between  cells  is  primarily  a  function  of  radial 
position.  Tube  powers  were  also  shown  to  have  a  radial  depen¬ 
dence  which  explains  why  02  concentration  profiles  within  power 
groupings  are  very  similar. 

3.4.3.  Contours  plots  and  cathode  oxidant  pathlines 

Contour  plots  of  bundle  temperature  and  surface  oxygen  con¬ 
centration  are  shown  in  Fig.  13.  In  both  figures  cathode  gases  enter 
through  the  dark  ring  at  the  bottom  of  the  bundle.  Cathode  gases 


flow  in  the  positive  z-direction  and  exit  the  bundle  through  concen¬ 
tric  circle  cutouts  in  the  outlet  tube-sheet.  Anode  gases  also  flow  in 
the  positive  z-direction  in  a  co-flow  configuration.  Symmetric  stack 
results  are  mirrored  to  show  half  of  the  bundle  in  Fig.  13.  Contours 
of  tube  temperature  are  shown  on  the  left  quarter  with  contours 
of  surface  oxygen  molar  concentration  shown  on  the  right  quar¬ 
ter.  Temperature  contours  illustrate  the  radiation  shielding  effect 
outer  periphery  tubes  have  on  inner  periphery  tubes.  The  inner 
periphery  ring  of  tubes  is  relatively  hot  with  a  uniform  tempera¬ 
ture  distribution.  Cold  zones  are  seen  at  the  outer  periphery  tubes 
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Fig.  13.  Contour  plots  of  surface  oxygen  molar  concentrations  on  the  right  quarter  and  tube  temperatures  on  the  left  quarter,  (a)  Looking  down  onto  interior  aspect  of  bundle, 
(b)  Looking  down  onto  exterior  aspect  of  bundle. 


Fig.  14.  Pathlines  with  cathode  colored  by  temperature.  Tube  bundle  and  central  fuel  preheat  tube  walls  shown  in  black,  (a)  Looking  down  onto  interior  aspect  of  bundle, 
(b)  Looking  down  onto  exterior  aspect  of  bundle. 


caused  by  their  close  proximity  to  the  stack  can  wall  as  well  as  by 
low  cathode  mixing.  As  each  tube  consumes  the  same  amount  of 
02  (constant  cell  current),  low  02  concentrations  at  outer  periph¬ 
ery  tubes  are  the  result  of  low  convective  mixing  not  02  diffusion 
at  the  cathode.  A  further  indicator  of  ineffective  cathode  mixing 
is  demonstrated  in  the  cathode  particle  pathlines  colored  by  tem¬ 
perature  in  Fig.  14.  To  reach  the  outer  periphery  tube  surfaces,  the 
oxidant  needs  to  turn  180°  after  entering  the  cathode.  Pathlines 
indicate  the  majority  of  oxidant  flows  interior  to  the  bundle  unable 
to  overcome  initial  inward  radial  momentum  and  reach  the  outer 
periphery  tubes. 

4.  Conclusion 

A  powerful  tubular  SOFC  system  design  and  simulation  tool  has 
been  developed.  A  detailed  stack  model  coupling  CFD  to  a  1-D 
electrochemical  tubular  cell  model  allows  one  to  see  how  stack 
geometries  affect  variations  in  tube  performance.  By  integrating 
a  recuperator,  tail-gas  combustor,  and  catalytic  partial  oxidation 
reformer  to  the  detailed  stack  model,  thermal  interactions  between 
BoP  and  the  SOFC  stack  are  captured  in  the  system  level  model. 


Model  results  can  point  fuel  cell  developers  to  more  effective  sys¬ 
tem  architectures  and  optimal  operating  conditions  with  the  goal  of 
increasing  system  efficiencies  by  optimizing  the  thermal  coupling 
between  BoP  and  the  SOFC  stack. 

The  model  capabilities  were  explored  through  simulation  of  a 
highly  integrated  tubular  SOFC  system  for  small  (~1  kW)  mobile 
applications.  The  simulation  predicts  a  fuel  conversion  efficiency 
of  21%,  where  75%  of  the  input  fuel  energy  is  convected  away  with 
the  exhaust  stream,  and  the  remaining  3%  lost  through  heat  transfer 
to  the  environment.  With  75%  of  energy  lost  in  the  exhaust  stream, 
model  results  highlight  the  ineffective  use  of  a  co-flow  recuperator, 
and  a  counter-flow  recuperator  is  a  recommended  design  improve¬ 
ment.  With  ineffective  preheating  of  oxidant,  cathode  gases  are 
148  °C  colder  than  anode  gases  entering  the  stack  causing  a  large 
decrease  in  tube  temperatures  near  the  cathode  inlet.  Unlike  planar 
stacks,  convective  cooling  is  not  the  dominate  mechanism  for  heat 
transfer  within  the  stack;  therefore,  a  large  AT  is  not  required  of 
cathode  gases  entering  the  bundle. 

Simulation  results  point  to  radiation  heat  transfer  as  the  dom¬ 
inate  mechanism  of  stack  cooling  within  tubular  stacks.  Radiation 
accounts  for  62-93%  of  total  heat  rejection  from  the  external  tube 
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surface.  The  dominance  of  radiation  leads  to  a  strong  relationship 
between  the  power  output  of  a  tube  and  the  view  factor  from  the 
tube  to  the  relatively  cold  stack  can  wall.  The  bundle  under  study 
results  in  seven  tube  groupings  based  on  similar  view  factors  to 
the  stack  can.  These  seven  view  factor  groupings  correspond  to 
seven  power  groupings.  Tube  performance  groupings  yield  insight 
into  potential  reduced-order  modeling  strategies  of  such  geomet¬ 
ric  configurations.  One  such  strategy  would  be  to  simulate  stack 
performance  based  on  tube  groupings  rather  than  the  extrapola¬ 
tion  of  a  single  electrochemical  tube.  This  strategy  would  require 
a  unique  set  of  thermal  boundary  conditions  to  be  extracted  from 
the  detailed  CFD  analysis  for  each  tube  grouping. 

An  enormous  amount  of  information  for  detailed  stack  design 
is  also  available  from  this  modeling  tool.  Contour  plots  of  stack 
temperatures  reveal  relatively  hot  and  uniform  tubes  at  the  inner 
periphery  but  cold  zones  develop  at  the  outer  periphery  of  the 
bundle.  Cathode  mixing  is  also  seen  to  be  relatively  low  with 
low  oxygen  concentrations  at  the  outer  periphery,  low  power 
tubes.  Non-uniformities  within  the  stack  lead  to  power  disparities 
amongst  cells.  As  cell  power  varies  from  7.6  to  10.8  W,  the  current 
stack  design  leaves  room  for  improvement.  It  is  recommended  that 
tubular  stacks  be  configured  such  that  the  variation  in  view  fac¬ 
tors  from  cells  to  the  stack  surroundings  is  minimized.  While  the 
simulation  tool  was  implemented  for  small-scale  tubular  SOFCs,  it 
should  be  noted  that  the  modeling  approach  employed  is  applicable 
to  a  wide  range  of  SOFC  systems. 
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Appendix  A.  TGC  model 

The  TGC  is  a  porous  annular  disk  where  unspent  fuel  from  the 
stack  is  oxidized  with  cathode  exhaust  gases.  TGC  exhaust  then 
enters  the  CFD  domain  at  the  recuperator  hot  flow  inlet.  As  with 
the  CPOx  model,  it  is  essential  to  couple  the  TGC  model  to  the  sys¬ 
tem.  As  an  example,  if  the  tube  model  predicts  low  stack  power, 
an  increase  in  unspent  fuel  will  increase  the  TGC  exhaust  tem¬ 
perature  which  indirectly  increases  the  temperature  of  air  leaving 
the  recuperator,  i.e.  entering  the  cathode.  A  higher  cathode  inlet 
temperature  acts  to  increase  tube  temperatures  and  power  out¬ 
put;  therefore,  TGC  coupling  acts  to  regulate  stack  power.  The  TGC 
model  domain  extends  from  the  top  of  the  outlet  tube-sheet  to  the 
top  of  the  system  insulation  (see  Fig.  2)  and  also  includes  the  anode 
gas  flow  through  the  inactive  tube  lengths  within  the  outlet  tube- 
sheet,  as  discussed  in  Section  2.4.4.  A  quasi  1-D  thermal  resistance 
model  created  in  Python  is  applied  to  the  TGC  domain,  as  shown  in 
Fig.  A.l,  to  capture  thermal  and  fluid  interactions  associated  with 
the  TGC.  Complete  combustion  is  assumed  within  the  TGC  [  1 ,3  ].  The 
TGC  domain  consists  of  a  mixing  region  where  stack  cathode  and 
anode  exhaust  gases  mix,  a  combustor,  combustor  exhaust  cavity, 
fuel/air  preheat  tube  flow,  air  preheat  tube  flow,  and  a  separating 
plate  that  separates  the  mixing  region  from  the  TGC  exhaust  cavity, 
as  shown  in  Fig.  6. 

A.  1.  Thermal  resistance  model 

Convection  and  conduction  heat  transfer  are  modeled  in  the 
same  manner  as  in  the  CPOx  model  where  all  surfaces  are  at  lumped 
temperatures  within  the  TGC  control  volume  (see  Fig.  6  for  geom- 
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Fig.  A.l.  TGC  model  thermal  resistance  network. 


etry).  The  surface  temperature  of  the  TGC  is  taken  as  the  average 
of  the  inlet  and  outlet  gas  temperatures.  An  effective  heat  transfer 
coefficient  is  used  in  the  TGC  domain  with  the  same  value  as  used 
in  the  CPOx  model  (see  Eq.  (10)).  A  perfectly  mixed  gas  condition 
is  applied  to  the  mixing  region  and  the  TGC  exhaust  gas  cavity  as 
given  in  Eq.  (12). 

A.2.  System  air  flow  preheat  modeling 

The  TGC  domain  also  has  pipe  flows  used  to  preheat  system 
streams  of  fuel/air  and  air.  Pipe  flow  is  not  perfectly  mixed,  and 
the  conductive  resistance  of  the  solid  pipes  is  assumed  negligi¬ 
ble.  There  are  four  air  preheat  flow  tubes  within  the  TGC  domain 
(only  two  seen  in  cut  view  of  Fig.  6).  A  turbulent  (Re  ^  5000)  pipe 
flow  Nusselt  number  relation  was  used  in  creating  a  mass  flow  rate 
functional  dependence  for  the  heat  transfer  coefficient  within  air 
preheat  tubes. 

hair,preheat  =  59664- m°f7  (A.l) 

where  mair  is  the  mass  flow  rate  (kg  s-1 )  of  air  per  air  preheat  tube. 

The  air  flow  tubes  pass  through  four  regions  of  the  TGC  domain, 
namely  the  insulation,  TGC  exhaust  cavity,  separation  plate,  and 
mixing  region.  The  magnitude  of  heat  transfer  from  each  region 
to  preheat  air  flow  is  calculated  by  combining  energy  balances 
and  rate  equations  applied  at  each  region  the  air  flow  tubes  pass 
through.  The  rate  of  heat  transfer  in  each  region  is  defined  as: 

Qair,i  —  ^air,preheat^air,i  [^avg,i  —  ^air,avg,i]  (A.2) 

where  Qair ,  is  the  thermal  energy  transferred  into  the  air  flow  trav¬ 
eling  through  tube  region  i  having  surface  area  Aair .  f  and  Tavg<i  is 
an  average  temperature  representing  the  air  flow  tube  in  region  i. 
When  the  tube  borders  a  solid  region,  Tavgi  is  the  average  between 
the  inner  and  outer  solid  surface  temperatures.  When  the  tube 
borders  a  gas  cavity,  Tavgl  is  the  perfectly  mixed  gas  tempera¬ 
ture.  Air  flow  temperatures  are  calculated  at  the  periphery  of  every 
region  that  the  flow  tube  intersects;  therefore,  Tair  avg  /  is  the  average 
between  the  inlet,  Tair  in  i,  and  outlet,  Tair  out ;,  air  temperature  within 
region  i.  Energy  balances  are  applied  to  air  flow  in  each  region  i  as: 

Qair,i  =  ^airQ^air  (^air,out,i  -  ^air,in,i)  (A.3) 
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where  Cpair  is  the  specific  heat  representing  air  within  the  preheat 
tube  and  is  evaluating  at  the  average  temperature  of  air  into  and 
out  of  the  entire  TGC  domain. 


is  the  specific  heat  of  the  mixture  calculated  at  the  average  outlet 
and  boiling  temperatures. 

The  rate  of  Qf/a,gas  is  calculated  with: 


A3.  System  fuel/air  flow  preheat  modeling 


Fuel/air  preheating  occurs  in  the  centrally  located  fuel  inlet 
flow  tube.  In  the  physical  system,  fuel  is  injected  (mixing  with 
air)  into  the  preheat  tube  by  an  atomizing  spray  nozzle  where 
the  nozzle  tip  is  located  at  the  outer  periphery  of  the  insulation. 
A  physical  spray  nozzle  is  not  modeled  within  the  TGC,  but  to  sim¬ 
ulate  this  boundary  condition,  a  mixture  of  liquid  fuel  and  gaseous 
air  enter  the  preheat  tube  at  a  common  temperature.  Within  the 
tube,  atomized  fuel  droplets  are  rapidly  vaporized.  Heating  of  the 
vaporizing  fuel  and  gaseous  air  mixture  is  modeled  by  considering 
the  air  and  liquid  fuel  components  separately.  An  energy  balance 
coupled  to  a  rate  equation  calculates  the  tube  area  where  vapor¬ 
ization  occurs.  Qf  ia  yap  is  the  total  energy  required  to  first  sensibly 
heat  liquid  fuel  and  gaseous  air  to  the  boiling  temperature  of  the 
fuel  and  the  latent  energy  required  to  completely  vaporize  the 
fuel. 


Qf/a,vap  —  ^Q^^boil  ^//a,in)ajr 

+  m  \hfg  +  Cpiiq  (Tboil  -  ^f  la, in)  ]  fuel  (A.4) 

where  Tboii  is  the  boiling  temperature  of  the  fuel,  T^a  in  is  the 
fuel/air  temperature  into  the  TGC  domain,  and  hfg  (3591<Jkg-1  for 
n-hexadecane)  is  the  latent  heat  of  vaporization.  The  mass  flow  rate 
of  air  and  fuel  entering  the  preheat  tube  are  mair  and  mfuel,  respec¬ 
tively.  Cpiiq  fuel  is  the  specific  heat  of  liquid  fuel  used  in  calculating 
the  sensible  heating  of  liquid  fuel.  Cpair  is  calculated  at  the  average 
of  the  inlet  and  boiling  temperatures. 

The  rate  at  which  Qf/a,vap  is  transferred  is  calculated  as: 


Qf/a, vap  —  ^boil^//a,vap 


/a,  tube  ' 


lf/a,  in 


+  T, 


boil 


(A.5) 


where  A^a  vap  is  the  surface  area  of  the  fuel/air  preheat  tube  from 
the  inlet  to  the  location  of  complete  fuel  vaporization.  T//a itube  is 
an  area  averaged  temperature  of  the  entire  fuel/air  preheat  tube. 
A  high  convective  heat  transfer  coefficient,  hboii  =  2000  Wm-2  K-1 , 
is  used  in  the  preheat  tube  as  a  simplified  means  to  simulate  fuel 
vaporization.  Rather  than  performing  a  detailed  analysis  involving 
non-dimensional  groups,  hboil  is  estimated  from  boiling  curves  [17]. 
This  rough  estimate  is  appropriate  because  system  level  predictions 
are  not  sensitive  to  hboil.  For  example,  an  hboil  =  1000  Wm-2  K_1 
only  decreases  the  fuel/air  temperature  leaving  the  TGC  domain  by 
0.7  °C  while  increasing  hboil  to  3000  Wrrr2  K_1  only  increases  the 
fuel/air  outlet  temperature  by  0.3  °C.  The  fuel/air  outlet  tempera¬ 
ture  is  not  sensitive  to  ftboil  because  of  a  low  convective  coefficient 
once  fuel  is  vaporized  as  discussed  in  the  following. 

After  the  fuel  has  completely  vaporized,  another  set  of  cou¬ 
pled  energy  balance  and  rate  equations  are  applied  to  determine 
the  extent  of  fuel/air  preheating  before  entering  the  CFD  domain. 
In  this  flow  region,  a  lower  convective  heat  transfer  coefficient, 
fy/a.pre  =  9.6  Wrrr2  K-1,  is  used.  hf/apre  is  based  on  a  laminar, 
Re  « 1 000,  pipe  flow  Nusselt  number  relationship  with  the  assump¬ 
tion  of  fully  developed  flow.  The  amount  of  sensible  heating  to  the 
fuel/air  gas  mixture,  Qf/a,gas>  is  calculated  as: 

Of  la, gas  —  thfiaCpf/a  (jf  fa,  out  —  ^boil)  (A.6) 

where  7}ya  out  is  the  temperature  of  fuel/air  gas  leaving  the  fuel/air 
preheat  tube.  The  combined  air  and  fuel  flow  rate  is  m^/a  and  Cpy/a 
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a,gas  —  tyva,gasA//a,gas 


J//a, tube  ' 


^boil  +  Tfi<j,  out 


(A.7) 


where  Af/a  gas  is  the  surface  area  of  the  fuel/air  preheat  tube  from 
complete  fuel  vaporization  to  the  outlet  of  the  fuel/air  tube.  It 
should  be  noted  that  the  area  averaged  temperature  of  the  entire 
fuel/air  preheat  tube  is  used  in  calculating  the  driving  potential 
in  Eqs.  (A.5)  and  (A.7).  The  reasoning  being  that  the  location  of 
complete  fuel  vaporization  is  dependent  on  system  level  operat¬ 
ing  parameters  and  TGC  model  dimensions.  With  the  location  of 
complete  vaporization  unknown,  it  is  difficult  to  calculate  an  area- 
averaged  fuel/air  preheat  tube  in  both  the  vaporizing  and  gaseous 
fuel  sections. 


A.4.  Model  integration  into  system  model 

The  thermodynamic  state  of  all  flow  inlets  entering  the  TGC 
model  which  include  the  cathode  and  anode  exhaust  are  extracted 
from  Fluent  and  the  tube  model  via  the  UDF.  The  state  of  preheated 
fuel/air  and  preheated  air  along  with  the  state  of  TGC  exhaust  leav¬ 
ing  the  TGC  domain  are  sent  as  inlet  conditions  to  the  CFD  model. 

As  in  the  CPOx  model,  there  are  also  solid  interfaces  at  the  TGC 
boundary  between  the  CFD  and  tube  models.  Adiabatic  boundary 
conditions  are  applied  at  the  interfaces  of  the  system  insulation, 
fuel/air  preheat  tube,  and  air  preheat  tube.  The  bottom  of  the  mix¬ 
ing  region  as  well  as  anode  gas  channels  in  inactive  tube  sections 
are  bound  by  surfaces  within  the  CFD  model.  At  these  CFD  sur¬ 
faces  a  convective  thermal  boundary  condition  is  applied  with 
h  =  1 00  W  rrr2  K-1  and  a  free  stream  temperature  equal  to  the  mix¬ 
ing  gas  temperature  as  calculated  in  the  TGC  thermal  resistance 
model.  The  total  heat  transfer  at  the  interface  calculated  by  the 
CFD  model,  Qcfdjgc,  is  added  to  the  energy  balance  at  the  mixing 
gas  node  within  the  TGC  thermal  model. 

References 

[1  ]  M.  Sorrentino,  C.  Pianese,  J.  Fuel  Cell  Sci.  Technol.  6  (2009). 

[2]  P.  Lisbona,  A.  Corradetti,  R.  Bove,  P.  Lunghi,  Electrochim.  Acta  53  (2007) 
1920-1930. 

[3]  N.  Lu,  Q.  Li,  X.  Sun,  M.A.  Khaleel,  J.  Power  Sources  161  (2006)  938-948. 

[4]  S.H.  Chan,  O.L.  Ding,  Int.  J.  Hydrogen  Energy  30  (2005)  167-179. 

[5]  L.  Petruzzi,  S.  Cocchi,  F.  Fineschi,  J.  Power  Sources  118  (2003) 
96-107. 

[6]  K.J.  Kattke,  R.J.  Braun,  J.  Fuel  Cell  Sci.  Technol.  8  (2)  (2011), 
doi:  1 0.1 1 1 5/1 .4002233. 

[7]  J.  Poshusta,  et  al.,  Solid  Oxide  Fuel  Cell  Systems  with  Hot  Zones  and  Two-Stage 
Tail  Gas  Combustors,  Protonex  Technology  Corporation,  US  Patent  Application 
12/006,668  July  9,  2009. 

[8]  A.M.  Colclasure,  B.M.  Sanandaji,  T.L.  Vincent,  R.J.  Kee,  J.  Power  Sources  196 
(2011)196-207. 

[9]  D.G.  Goodwin  An  open-source,  extensible  software  suite  for  CVD  process  sim¬ 
ulation,  in:  M.  Allendorf,  F.  Maury,  F.  Teyssandier,  (Eds).,  Chemical  Vapor 
Deposition  XVI  and  EUROCVD  14,  volume  PV  2003-08,  pages  155-162.  Elec¬ 
trochemical  Society,  2003.  see  also  http://code.google.eom/p/cantera/. 

[10]  Special  Metals,  INCONEL  600  Data  Sheet,  product  data  sheet,  Table  3,  p.  2, 
September  2008. 

[11]  Microtherm  International  Ltd.,  Microtherm  Insulation  Product  and  Perfor¬ 
mance  Data,  product  brochure,  Table  3,  p.  5,  June  2001. 

[12]  F.  Calise,  M.D.  d’Accadia,  G.  Restuccia,  Int.  J.  Hydrogen  Energy  32  (2007) 
4575-4590. 

[13]  J.  Jia,  A.  Abudula,  L.  Wei,  R.  Jiang,  S.  Shen,  J.  Power  Sources  171  (2007)  696-705. 

[14]  G.  van  Rossum,  Python  (Version  2.6)  [Computer  Software].  Available  from: 
http://www.python.org. 

[15]  J.  Krummenacher,  K.  West,  L.  Schmidt,  J.  Catal.  215  (2003)  332-343. 

[16]  K.  Hohn,  T.  DuBois.J.  Power  Sources  183  (2008)  295-302. 

[17]  F.  Incropera,  D.  DeWitt,  T.  Bergman,  F.  Lavine,  Fundamentals  of  Heat  and  Mass 
Transfer,  Fig.  10.4,  6th  ed.,  Wiley,  New  York,  2007. 


